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We present a general formula of the orbital magnetization of disordered systems based on the 
Keldysh Green's function theory in the gauge-covariant Wigner space. In our approach, the gauge 
invariance of physical quantities is ensured from the very beginning, and the vertex corrections are 
easily included. Our formula applies not only for insulators but also for metallic systems where 
the quasiparicle behavior is usually strongly modified by the disorder scattering. In the absence 
of disorders, our formula recovers the previous results obtained from the semiclassical theory and 
the perturbation theory. As an application, we calculate the orbital magnetization of a weakly 
disordered two-dimensional electron gas with Rashba spin-orbit coupling. We find that for the short 
range disorder scattering, its major effect is to the shifting of the distribution of orbital magnetization 
corresponding to the quasiparticle energy renormalization. 

PACS numbers: 75.10.Lp, 73.20.Hb 



I. INTRODUCTION 

Magnetization is one of the most important and in- 
triguing material properties. An adequate account of 
magnetization should not only include the contribution 
from the spin polarization of electrons, but also the con- 
tribution from the orbital motion of electrons. In crys- 
tals, due to the reduced spatial symmetry, the orbital 
contribution to the magnetization is usually quenched. 
However, in certain materials with topologically nontriv- 
ial band structures, large contributions can arise from the 
effective reciprocal space monopoles near the band anti- 
crossings. Several different methods have been employed 
to study the orbital magnetization (OM) in crystals 1 - - — . 
One major difficulty in the calculation is posed by the 
evaluation of the operator 1,2 f x j, because the the po- 
sition f is ill-defined in the Bloch representation. This 
difficulty can be avoided in a semiclassical picture or be 
circumvented by a transformation to the Wannier repre- 
sentation. Xiao et alA& presented a general formula for 
OM for metal and insulator, derived from a semiclassical 
formalism with the Berry phase corrections. Thonhauser 
et al&Z- derived an expression of the OM for periodic in- 
sulators using the Wannier representation. From the ele- 
mentary thermodynamics, Shi et alP obtained a formula 
for the OM in a periodic system using the standard per- 
turbation theory. Their result can in principle take into 
account the electron-electron interaction effects. A com- 
putation of the OM for periodic systems with density- 
functional theory was carried out by Ceresoli et alJ^. 

Previous studies are mainly concerned with clean sys- 
tems. However, real crystals are never perfect, disor- 
ders such as defects, impurities, phonons etc. constantly 
break the translational symmetry and lead to scattering 
events. The effect of disorder scattering on the OM has 
not been carefully studied so far. On one hand, the OM 



is a thermodynamic quantities, hence it is expected to 
be less susceptible to disorder scattering. On the other 
hand, the appearance of current operator j in the defi- 
nition suggests behaviors similar to transport quantities 
which might be strongly affected by the disorder scat- 
tering. Therefore, it is important and desirable to have 
a good understanding of the role played by the disorder 
scattering in the OM. 

In this paper, we present a general formula of the OM 
in disordered systems based on the Keldysh Green's func- 
tion theory in the gauge-covariant Wigner spaced— . 
This approach was developed as a systematic approach 
to the nonequilibrium electron dynamics under external 
fields. Our formula derived from this approach shares 
the advantage of being able to capture the disorder ef- 
fects in a systematic way and ensure the gauge invariance 
property from the very beginning. We show that in the 
clean limit, our formula reduces to the previous results 
obtained from other approaches. As an application, we 
study the OM in a disordered two-dimensional electron 
gas with the Rashba spin-orbit coupling. We find that 
the OM is robust against short range disorders. The main 
effect of the scattering by short range disorders is a rigid 
shift of the distribution of OM in energy. 

The structure of this paper is organized as follows. In 
Sec. Ill Al we outline the Keldysh Green's function for- 
malism which is employed for our derivation. Our general 
formula of OM is presented in Sec. Ill Bl In Sec. IIIII we 
apply the formula to study the OM of a two-dimensional 
disordered electron gas with the Rashba spin-orbit cou- 
pling. Summary and conclusion are given in Sec. IIVI 
Some details of the calculation are provided in the ap- 
pendices. 
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II. ORBITAL MAGNETIZATION OF 
DISORDERED SYSTEMS 

A. Keldysh Green's function formalism 

We employ the Keldysh Green's function formalism 
in the Wigner representation^, which has recently been 
used to study the current response of multi-band systems 
under an electric fiel d 19 ' 20 . In the Wigner representation, 
Green's functions and the self-energies are expressed as 
functions of the center-of-mass coordinates (T,X), the 
energy e and the mehanic momentum p. The energy 
and the mechanic momentum are the Fourier transforms 
of the relative time and space coordinates respectively. 

The Dyson equations in the presence of external elec- 
tromagnetic fields can be written as 



*<L-Eo(p)-Ue) *G(e,p) = I 



G(e,p)* eL-Ho(p)-Ue) 



(la) 
(lb) 



Each quantity with an underline in the above equations 
is a matrix in Keldysh space. Specifically, we have 
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where & RA< ^ are the (retarded, advanced, lesser) 
Green functions, and s'- R >- A > < ) are the corresponding self- 
energies, Hq is the Hamiltonian in the absence of external 
electromagnetic fields, <t° is the identity matrix. The * 
operator in Eq.([T]) is defined as 



exp 



(4) 



with the differential operators d and ~~§ operating on the 
left-hand and the right-hand sides respectively, q = — e 
is the electron charge, and = d Xli A v {X)-d Xv A»{X) 
is the electromagnetic field tensor, [i and v label the 
four dimensional space-time components and the Ein- 
stein summation convention is assumed. It should be 
noted that the energy e and the mechanic momentum p 
include the electromagnetic potentials A^(X), both are 
gauge invariant quantities. The * operator in Eq.([T]) only 
involves the physical fields, so it is also gauge invariant. 
In this formalism the gauge invariance is respected from 
the very beginning and easily maintained during the per- 
turbative expansion, which is an important advantage^. 

Here we consider the situation with a uniform weak 
magnetic field along the z-direction, i.e. B = (0,0, B). 
Then the various quantities can be expanded in terms of 
B. In particular, Green's functions and the self-energies 



can be expressed as 

G a {e, p) = Gff(e, p) + eHBG%(e, p) + 0{B 2 ), (5) 
t a (e)^±^s) + ehBt%(e) + 0(B 2 ), (6) 

with a — R, A, < for the retarded, advanced and lesser 
components respectively. Here functions with the sub- 
script are of zeroth order in the external magnetic field 
strength (note that they include scattering effects). We 
have 



G 



R(A) 



(e»P) 



e-H Q (p)-£ R(A \e) 



(7) 



G<(e,p) = G A (e,p)-G R (e,p) /(e), (8) 



where f(e) is the Fermi distribution. The functions with 
subscript B are the linear response coefficient to the ex- 
ternal field. They can be solved from the Dyson equation. 
It is usually convenient to decompose the lesser compo- 
nent Gg and E^ (which are related to particle distribu- 
tion) into two parts, with one part from the Fermi surface 
and the other part from the Fermi sea^, 

G%{e,p)= £<>,p)fc/(e}+<5< (e,p)/(e) (9) 



E<( £ ) = S< (e)d e f(e) + t< (s)f(e) 



(10) 



From the Dyson equation (kept to the linear order in B), 
it is straightforward to show that 



^BJ 



0, 



(11) 



i.e. there is no Fermi surface term in the linear order 
lesser component, and for the Fermi sea term we have 



U B,II 
n B.ll 



(e,p) = G A (e,p)-G R (s,p), 
(e) = ± A (e)-± R (s). 



(12) 
(13) 

The retarded and advanced Green's function Gg^ and 

self-energy E^^ are determined from the following self- 
consistent equations 



G 



*,R(A) _ 1 



G R(A) v x (d Py G^>) - (d Py GT>)v x G» 



yR(A) 



-G. 



R(A)f,R(A)gR(A) 











(14) 



where the velocity operator is defined as Vi = i 



Xi,H 



In this approach, the disorder effects arc captured by 



the self-energies E w and J^^ A> , which allows a sys- 



*< A > and tf A \ 
tematic perturbative treatment. In the weak disorder 
regime, the self-consistent T-matrix approximation pro- 
vides a good approximation scheme. In this approxima- 
tion, we have 



\R(A) 



(e) = n imp f R{A) (e), 



(15) 



and 



\R(A) 



(e)=n 



T 



R(A) 



00 



d 2 p 

(2irhy 



G$ A \e,p)T^'(e), 



R(A), 



(16) 
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where ni mp is the impurity concentration and the T- 
matrix is expressed as 



T, 



(e) = V in 



1 - 



(27rft) ; 



imp 



(17) 



with Mmp being the impurity potential. 

The equilibrium Green's functions G^ A ^ and the self 

energies can be obtained by solving Eqs. 0, (fT5|) 

and (jTTJ) self-consistently. Then the linear order coef- 



ficients and llg (A> can be solved from Eqs. (TTJ} 



and (|16p . Finally, we can obtain Gg n through Eq. (TTJ 
and the linear response of the system in the external mag- 
netic field can be completely determined. 

The lesser Green's function contains the information of 
particle distribution. In our case, both the external mag- 
netic field and the disorder scattering affect the quasi- 
particle distribution. Before we proceed, it is interesting 
to observe how the non-trivial band geometry (described 
by the Berry curvature) can be captured by the present 
Wigner space Green's function formalism. For a homo- 
geneous system, the electron density can be written as 



1 



ds 
2^ 



d 2 p 

{2irh)- 



■tr 



G<(e,p) 



(18) 



In the absence of the disorder scattering, the eigenstates 
are well-defined Bloch states grouped into energy bands. 
Using the theorem of residues, we can express the ground 
state electron density in the presence of a constant mag- 
netic field as (see Appendix [C]) 



E 



d 2 P 



(19) 



The summation is over all the occupied states, and 
Cl n (p) = «(V p u„p| x |V p M np ) is the Berry curvature of 
the Bloch state |n, p) = e* p ' x / h \u np ). It can be seen that 
the Fermi-sea volume is changed linearly by a magnetic 
field when the Berry curvature is nonzero. This effect was 
previous interpreted as the modification of phase space 
density of states^. 



B. Formula of orbital magnetization 

We start from the standard thermodynamic definition 
of the OM density at zero temperature^: 



M 



dK 
~dB 



(20) 



where K = E — /iN is the grand thermodynamic po- 
tential, B is a weak magnetic field. Since we are con- 
cerned with the orbital contribution, the small Zeeman 
coupling between the electron spin and external field will 



be dropped. The potential K can be expressed through 
the lesser Green's function, 



K 



(21) 



Using Eqs. (HOI), (ED) and ©, we find that the OM can 
be written as 



M = -ieh / — 
2?r 



d 2 p 



tr 



(27rft) 2 
G A (s,p)-G§(e 7 p 



fie)- 



(22) 



From this expression, we can see that the OM has con- 
tributions from the whole Fermi sea, with no separate 
Fermi surface contribution such as that for the transport 
quantities. 

In this formula, the impurity scattering effect comes 
in through two terms: the self-energy Y,q' A which modi- 
fies the ground state electronic structure and the vertex 
corrections associated with Yi^ A,< which represent an 
interplay between the magnetic field and the impurity 
scattering. We may separate out the terms containing 
E B ' ' and write the OM explicitly as 



where 
M 1 



eh 

T 



M = M L +M 



— f(\[ d2p 
2ir n£, J {2nhf 



II 



(23) 



Tr [eii(H - ^)G A (e, p)viG A (e, p)vjG A (e, p) 



(24) 



and 



M u = -ieh 



ds 



xTr 



2tt 

^)G a E^Gq 



d 2 p 

(2nh) 2 



(25) 



where with i,j 6 {x, y} is the 2D antisymmetric ten- 
sor, and the second term in the bracket in Eq. (|24[) means 
that the second term is the same as the first term except 
that all the Gg are replaced by Gq. Such a decomposi- 
tion scheme was also adopted in the study of anomalous 
Hall conductivity^, and in that context, the two parts 
are referred to as the intrinsic part and extrinsic part re- 
spectively. It should be noted that the intrinsic part M 1 
also has impurity scattering effects in it (see Eq. ([7]) and 
Eq. (|15p ). it is intrinsic in the sense that it only contains 
quantities that are of zeroth order in the external field. 
As for the extrinsic part M 11 , it is easy to see that it is 
already linear order in n- lm p (see Eq. (|16[) ). Therefore in 
the weak scattering regime, the extrinsic part is expected 
to be much smaller than the intrinsic part. 

The above formula is our main result. From this for- 
mula, we see that there is no separate Fermi surface con- 
tributions like those in the transport quantities, which is 
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consistent with OM being a thermodynamic equilibrium 
property. This formula applies for both insulators and 
metals. The quantities in this formula can be calculated 
from the Dyson equation according to our prescription 
described in the previous section. It can also be straight- 
forwardly implemented in the numerical calculation, ei- 
ther from effective models or from first principles. 

In the clean limit, we only have the intrinsic part. The 
general result reduces to (see Appendix [D] for the deriva- 
tion) 



M = ^2 fnp m n(p) - r(e np 

7ip 



/i)tt„(p) , (26) 



moment of the Bloch state \n, p) and 
x | V D M ra p) is the Berry curvature. The 



where m n {p) = {e/2h)i(V p u np \[e n (p)-H Q (p)] x |V p w np ) 
is the orbital 

fin(p) = i(V p U„p| A|, p 

first term in Eq. (1261) is a sum of the orbital magnetic mo- 
ments associated with each Bloch stat o 24 i 25 , and the sec- 
ond term is a Berry-phase correction to the OM. There- 
fore, the OM can be written as 



M = Mr 



(27) 



This clean limit result was previously derived from the 
standard perturbation theory of quantum mechanics by 
Shi et al. — and also from the semiclassical theory by Xiao 
et al. — . Now it is also reproduced as a special limiting 
case of our general formula. 



III. APPLICATION TO A TWO-DIMENSIONAL 
ELECTRON GAS WITH RASHBA SPIN-ORBIT 
COUPLING 

A. Model 

We apply our theory to study the model of a two- 
dimensional disordered electron gas with Rashba spin- 
orbit coupling. The Hamiltonian for the system reads 
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FIG. 1. (Color online) (a) Electronic band dispersions of our 
model given by Eq. (|29l) . (b) Orbital magnetization (M, solid 
red curve) of disordered free system and its two components 
Mm (dashed blue curve) and Mn (dash-dotted green curve) 
as functions of Fermi energy Ef- They are plotted in units 
of e/h. The parameters are chosen as 2ma 2 — 3.59 and Ao = 
0.1. 



where A = 1 , 2 labels the upper and lower band respec- 
tively. When the Rashba coupling energy scaled ma 2 
is larger than the Zeeman coupling strength the minima 
of the lower band occur at a finite wave vector and the 
dispersion assumes a Mexican hat shape (see Fig.[T](a)). 
When the Zeeman coupling dominates over the Rashba 
energy, the minimum of the lower band Ei occurs at the 
origin (see Fig. [2] (a)). 

Let's first consider the clean limit, in which case the 
orbital magnetic moment and the Berry curvature of each 
Bloch state can be calculated straightforwardly: 



H — H o + i?imp, 



(28a) 



2 

= ^-^ + a ( Px a y - Py a x ) - A a z , (28b) 
2m 



^imp — ^imp 



(28c) 



where (a x ,a v ,a z ) are the three Pauli matrices and cr° 
is the identity matrix, a is the strength of the spin-orbit 
coupling, and term — Aq& z is the spin splitting which can 
be introduced by the exchange coupling with a nearby 
ferromagnet or magnetic dopants. i?i mp is the disorder 
potential from the randomly distributed short range im- 
purities with strength Wi mp . The energy dispersion of the 
Hamiltonian Hq is given by 



Exip) 



P 

2m 



a 2 p 2 , 



(29) 



mi(p) = m 2 {p) 



fii(p) = -fi 2 (p) 



A a 2 



2h A§ 



2 ? ' 

a p 



A a 2 



2(A 2 +a 2 p 2 )i 



(30) 
(31) 



It is interesting to observe that for the same wavevec- 
tor the orbital moments of the two bands have the same 
magnitude and the same sign, while the Berry curvatures 
have the same magnitude but opposite signs. It should 
also be noted that both the orbital moment and the Berry 
curvature would vanish if either a or Ao vanishes. From 
Eq. (|26|) . we further see that the OM is nonzero only when 
both the spin-orbit coupling and the exchange coupling 
are present. 

Analytical expressions of the OM can be easily ob- 
tained for the clean limit using Eg. (1261) . For example, for 
the case with Ep > Aq, we have 
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FIG. 2. (Color online) (a) Electronic band dispersions of our 
model given by Eq. (|29[) . (b) Orbital magnetization (M, solid 
red curve) of disorder free system and its two components 
Mm (dashed blue curve) and Mn (dash-dotted green curve) 
as functions of Fermi energy Ef- They are plotted in units 
of e/h. The parameters are chosen as 2ma 2 — 0.08 and Ao = 
0.1. 




FIG. 3. (Color online) (a) Orbital magnetization (OM) as 
a function of Fermi energy Ef with different impurity con- 
centration riimp. (b) Density of OM with different impurity 
concentration ni mp . These quantities are plotted in units of 
e/h. The parameters are chosen as 2ma 2 — 3.59, Ao = 0.1, 
and uimp = 0.1. 
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a P p 2 



8Trmha 2 

where pp t 2 is the Fermi momenta of the two bands. 
B. Results 

Now we analyze the OM of the disordered 2D Rashba 
model in detail. The calculation procedure follows our 
discussion in Sections III Al and III Bl Since we have seen 
that both the spin-orbit coupling and the exchange cou- 
pling are essential ingredient for the OM, in the following 
we shall consider two different regimes of the model deter- 
mined by the competition between the Rashba spin-orbit 
coupling and the exchange coupling. For each regime, we 
first analyze the clean limit where the physical picture is 
more transparent, and then study the influence of disor- 
der scattering which is the focus in this paper. 

We first consider the regime where the Rashba coupling 
dominates over the exchange coupling, i.e. 2ma 2 ^> Ao. 
The typical band dispersion in this regime is shown in 
Fig. [1] (a) (with 2ma 2 = 3.59 and A = 0.1). In this 
regime, the bottom of the lower band occurs at a finite 
wavevector. The energy spectrum around the origin has 
an effective Dirac cone structure with a local gap 2Aq at 
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FIG. 4. (Color online) (a) Orbital magnetization (OM) as 
a function of Fermi energy Ef with different impurity con- 
centration nimp. (b) Density of OM with different impurity 
concentration ni mp . These quantities are plotted in units of 
e/h. The parameters are chosen as 2mce 2 — 0.08, Ao = 0.1, 
and uimp = 0.1. 
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p = 0. Both the orbital moment and the Berry curva- 
ture are concentrated near this band anticrossing point, 
as is evident from Eqs. (f30"l) and (|3Tj) . Fig. Q] (b) shows 
the OM for the clean limit. The orbital moment contri- 
bution A4 m and the Berry curvature contribution A4q 
are also plotted in Fig. Q] (b). We can see that as the 
Fermi energy Ep increases from the lower band bottom, 
M.Q, increases while M. m decreases. The increasing rate 
of Mn is higher than the decreasing rate of Mm, so the 
overall OM is increasing. The OM reaches its maximum 
when Ep — — A , which corresponds to the local band 
top around the origin in momentum space. As the Fermi 
energy sweeps across the local energy gap between — Ao 
and +Ao, the OM decreases approximately linearly with 
Ep. The linearity can be understood by noticing that 
from Eq.([2H])) the derivative of the OM with respect to 
Ep is just the momentum space integral of the Berry 
curvature. The Berry curvature distribution is concen- 
trated near the band anticrossing point, corresponding to 
the small region around the origin in the present model. 
When the Fermi energy is within the gap, the Berry cur- 
vature integral only has contribution from the lower band 
and is almost constant, therefore leading to the linear 
energy dependence of OM. This linear decrease of OM 
stops when the Fermi energy touches the bottom of the 
upper band at +Ao. Above the upper band bottom, Mq 
and M m almost cancel each other and the OM is vanish- 
ingly small. Throughout the spectrum, M m is positive 
while Mn is negative, corresponding to the paramag- 
netic and diamagnetic responses respectively. This has 
a clear explanation in the semiclassical picture: Mm is 
due to the self-rotation of the wavepacket which is para- 
magnetic, while Mn is from the center-of-mass motion 
of the wavepacket hence is diamagnetio^. 

When the exchange coupling dominates over the 
Rashba energy, The minimum of the lower band occurs at 
the origin. Compared with the previous case, there is no 
local gap at p = 0. The typical band dispersion is shown 
in Fig. [2] (a) (with A = 0.1 and take 2ma 2 = 0.08). 
The overall shape of the OM is similar to that for the 
first case. Its distribution over spectrum is mainly below 
the upper band bottom. However, due to the absence of 
the local gap, the kink point at — Ao in Fig.[T](a) merges 
with the lower band bottom. Moreover, the two contri- 
butions Mq and M m strongly cancel each other and the 
resulting OM is much smaller. 

Now let's consider the effects of disorder scattering on 
the OM in our model. When the disorder scattering is 
turned on, the translational invariance is broken. We can 
no longer define quantities such as Mq, and Mm- Their 
effects are merged into the sophisticated expression in 
Eq.([23]). Fig. (a) and Fig. H (a) show the OM versus 
Ep for the two regimes we discussed above. The different 
curves in each figure correspond to different impurity con- 
centrations 7ii, np . Compared with the clean limit where 
"■imp = 0, we see that the shape of the OM curve is al- 
most unchanged but mainly its position is shifted by the 
scattering. This behavior is more obvious when we look 
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FIG. 5. (Color online) Orbital magnetization (OM) as func- 
tions of Fermi energy Ef and the impurity concentration rii mp 
in units of e/h. The parameters are chosen as Ao = 0.1, and 
Wimp = 0.1, except the Rashba energy: (a) 2ma 2 = 3.59, (b) 
2ma 2 = 0.08. 

at the density of OM shown in Fig. (b) and Fig. H (b). 
For the clean limit, we see that the major contribution to 
the OM is from the states at the band bottom and at the 
local band edge. The effect of disorder scattering here 
is to shift the the density of OM distribution in energy. 
Such a shift can be understood by noticing that the OM 
only has the Fermi sea contribution. The main effect of 
scattering in Eq. (l23l) is the shift of energy arising from 
the real part of the self-energy correctional. For the short 
range disorder model, the disorder potential is a constant 
in momentum space, hence the sclf-cncrgy is independent 
of the state, which results in a rigid energy shift for all the 
states. For a general disorder potential, the energy shift 
would be generally different for different states therefore 
the distribution of OM would be distorted. The effects of 
finite range disorders are currently under investigation. 

To leading order, the shift should be linear in the disor- 
der density n mip . In Fig. [5] we plot the OM as a function 
of Ep and 7ii mp . The linear dependence of the energy 
shift in rtj mp is clearly observed. Apart from the energy 
shift, the scattering induced state broadening is mani- 
fested as the smoothing of the peaks of the density of 
OM, which can be clearly observed in Fig. [3] (b) and 
Fig. 2] (b). The peaks of OM are only slightly decreased 
by the scattering. This means that the OM carried by 
the electronic states are quite robust against scattering. 

IV. CONCLUSIONS 

In summary, we have derived a formula of the the 
OM of disordered electron systems based on the Keldysh 
Green's function theory. This approach was developed 
as a systematic approach to the nonequilibrium electron 
dynamics under external fields. In the formula, OM is 
expressed in terms of the Green's functions and self- 
energies, which can be solved from the Dyson equations, 
and systematic approximation schemes to the disorder 
effects can be employed. We find that there is no Fermi 
surface contribution like in the case of the current re- 
sponse. Our formula applies not only for insulators but 
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also for metallic systems, where the quasiparicle behav- 
ior is usually strongly modified by the disorder scatter- 
ing. It can also be straightforwardly implemented in the 
numerical calculation. In the clean limit, our formula 
reduces to the previous result obtained from other ap- 
proaches. As an application, we calculate the OM of 
a weakly disordered two dimensional electron gas with 
Rashba spin-orbit coupling. The result shows that in the 
simplest white noise short range disorder model, the OM 
is robust against weak scattering and the main effect of 
scattering is a rigid shift of the distribution of OM in en- 
ergy, which can be attributed to the real part of the self 
energy. 
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Appendix A: Self-consistent equation for T,q and 
explicit forms of Go 

The Green functions and self-energies in the absence 
of the external fields are obtained from the coupled self- 
consistent equations (J7J, (|T5j) and p7|) . In our model, 
a direct analytical integration in e shows that 



and 



/,>!! ^imp^imp (l ^imp^o (^)) 

(l-u imp g$°(e)) 2 -u 2 nip g^(e) 2 



, (Al) 



^ z{£) = (i - u imp9 f(e)y- ^g^W ' (A2) 

££*( £ ) = E*»( e ) = 0, (A3) 



50 (£) 47TS 2 ^' li G$(e,0,a) 



ma 2 g*(e), (A4) 
(A5) 



+ H- Z™(e) + ma 2 + *R R (e)) ]^,(A6) 



G*(e,p,±) = (e-j^m+M-ETO 



T V«V + (-Ao + ^ 2 (e)) 2 ) ,(A7) 
R R (e) = ((ma 2 ) 2 + 2ma 2 (e + /x - E£°(e)) 

+ (-A + Ef (e)) 2 Y, 

(A8) 

where A is the cut-off in momentum integration, and 

G R (s) = G R0 (e)a° + £ G Rl (s)a l , (A9) 

I— x,y,z 

with 

G R0 (e, p) = (e - p 2 /2m + M - £*°( e )) G R (e,p), 

(A10) 

G™(e,p) = (-ae ijzPj + £ fa (-A + E**(e))) G R (e,p), 

(All) 

G R (e,p) = (e - p 2 /2m + p - E™(s)) 2 

+ aV + (-A Q + Ef ( £ )) 2 , (A12) 

and eiji is the anti-symmetric tensor, (i, j, I, • • • ) label 
the Cartesian components. The same results have been 
obtained in Ref. I2CJ. 

For each e, the self-energy can be calculated by iter- 
ations which can be performed until the the prescribed 
accuracy is reached. 



Appendix B: Self-consistent equation for Gg and 
and their explicit forms 

The equations for solving the first order corrections Gg 
and E§ are presented here. Using Eqs. (fT4"|) and (TTB|). the 
retarded Green's function Gg can be rewritten as 



G R ( £ ) = G R0 (e)a + G R (e)-a, 



(Bl) 



with 



G R0 (e, p) = (G R0 (s, p) 2 + G R (s, p) 2 )E™( £ ) + 2G™(e)G«(e, p) • £* (e) 

+ G«( £ ,p)(d Pi tf (p) x G R (e,p)) • (^,ff (p)), (B2a) 
G R (e,p) = G «(e,p)Ei(e) - G *(s, p)G R0 (e, p)(d P Mp)) x (9 Py # (p)) 

-G*(^ff Q °(p))G«( e ,p) x (d Py H (p)) + G R (d Py H° (p))G R (e,p) x &.flo(p)) 

+ 2G«( £ ,p) (G «°(e,p)Ei°( £ ) + G«(e,p) • £|( E )) , (B2b) 
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and the inner product of two vectors are denned as 



A-B= J2 ^ Bl - ( B3 ) 



From Eqs. (fT6|) and (|17[) . we write the self-energy as 



with 



Sg(e) = S™( £ )a + (B4) 



Ef( £ ) = n lmpU ? mp ((1 ~ "imp3o fl0 (e)) 2 - uf mp g^(s)Y 2 

x [((1 - u imp5o flo ( £ )) 2 + u\ mp g**{sf) gf(e) + 2(1 - u imp g™(e))u imp g^(e)g« z (s)] , (B5a) 

= n impW f mp ((1 - u imp5o «°(e)) 2 - uf mp g^(ery 2 

x [((1 - u imp g™(e)f + uf mp g^(e) 2 ) g**(e) + 2(1 - u imp g™(e))u imp g** ( S )g™(e)] , (B5b) 
E*(e) = n impU f mp ((1 - Uimp g fl0 ( £ )) 2 - uLpSf*^) 2 )" 1 (B5c) 



9o,b& = / 7^2 G o% (B6) 



and we have 

, (27Tft)2 G °> B ' 

where a 6 {0, x, y, z}. The zeroth order components Gq 01 are computed as in appendix lAl and are used as input for 
the above equations. 

Appendix C: The particle density 

Here, we present the derivation of Eq. (|19[) . In the absence of disorder scattering, 

G«- A (e,p) = {e-H a (p)±i0 + }-\ (CI) 
At zero temperature, plugging Eq. (jCll) into Eqs. (fT2t and (|T8l) . we can obtain 

fde f d 2 p fy* 1 y, 1 1 

" e " 7 7T 7 (27rft) 2 I ^ e - e„ P + i0+ + & ^ (e - e„ p + i0+) 2 e - e mp + i0+ 

x ^[{unp\i> x (p)\u mp }{u mp \v y (p)\u np }}} . (C2) 

u np are the eigenfunctions of the unperturbed Hamiltonian and e np the eigenvalues. The integral over e contains 
simple and double poles. Using the residue theorem 2 ^, we obtain 

f d 2 

n e = I ^ ^2 ^{ 1 + 2iehB ^^[{ u n P \vx{p)\u mp }(u mp \v y (p)\u np }}}, (C3) 



II .()(-<- 



where occ denotes summing over occupied states. Further simplification can be made by using the Sternheimer 
equation 

Vj(p)\u„ p ) = (enp - e n 'p)\-^p-) + -^-\u np ), (C4) 
and we finally arrive at the equation 

d2p - h + £i 
h 



"•=E/pi i + s B-nn(p) ■ (C5) 
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Appendix D: Orbital magnetization in the clean limit 



The derivations of Eq.(f27| for the OM in the clean limit are present below. When the relaxation rate vanishes, 
substituting Eq. (jCll) into Eq. (|2~4l . we can write Eq. as 

M = eh 1 7^/(e) j j^^^( e np- ^)^{u„ P \^p)\u mp ){u mp \vy{p)\u np )] 

" V / nr. m 



1 



. (e - £n P + i0 + ) 2 e - e mp + i0+ (e - e„ p + i0+) 2 e - e mp + i0+ 
Using the residue theorem, we find that 

d 2 P „,/(w)-/( e np) . /'(e« P ) 



(Dl) 



M = -eh 



(2ttH) 2 (e mp - e„ p ) 



x y 



(p)\u np )} (D2) 



where f' np = df(e np )/de np . With the help of the Sternheimer equation Eq. (|C4j) . we obtain 



M = -eh 
2 



rf 2 p 

(27r/l)< 



E 



-p - PA-dpl^P - #o(p)] x |^E)/; p - (^gP|[ £np + # (p) - x l^gP-^p 



The above result can be written as 

M = ^2 { TO «(P)/»P + ( £ ™P ~ M) m "(P)/np - |( e rip - Ai)^np(p)} 



(D3) 



(D4) 



where m„(p) = (e/2?i)i(VpU, i p| [e„(p) — _ff (p)] X |V p u np ) is the orbital moment of state n, p and fi n (p) = 
i{V p u np \ x |VpU„ p ) is the Berry curvature. At zero temperature, /' becomes a <5-function of (e„ p — /x), therefore we 
have in this case 



M = ^2 [ m «(P)/«P ~ ^( e «P ~ M) fi n(p) 

np 



(D5) 
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